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Python Classes for Numerical Solution of PDE’s 

Asif Mushtaq, Member, lAENG, Trend Kvamsdal, Kare Olaussen, Member, lAENG, 


Abstract —We announce some Python classes for numerical 
solution of partial differential equations, or boundary value 
problems of ordinary differential equations. These classes are 
built on routines in niunpy and scipy. sparse . linalg (or 
scipy.linalg for smaller problems). 

Index Terms —Boundary value problems, partial differential 
equations, sparse scipy routines. 

I. Introduction 

T he Python computer language has gained increasing 
popularity in recent years. For good reasons; It is fast 
and easy to code and use for small “prototyping” tasks, since 
there is no need for explicit declaration of variables or a 
separate compilation cycle. It is freely available for most 
computer platforms, and comes with a huge repository of 
packages covering a large area of applications. Python also 
have features which facilitates development and encourages 
documentation of large well-structured program systems. 

Obviously, as an interpreted language native Python is not 
suitable for performing extended numerical computations. 
But very often the code for such computations reduces to 
calls to precompiled library routines. The numpy m and 
scipy lEl, llll packages make a large number of such 
routines directly available from Python. These packages are 
freely available for most operating systems, including Linux, 
OSX, and MSWindows. 

We here describe a process of making some of these 
routines even simpler to use for a field of applications, the 
numerical solution of partial differential equations discretized 
on a rectangular grid (or a subdomain of such a grid). As a 
simple reference problem one may consider the solution of 
the wave equation in the frequency domain, 

{-A + uj'^) (p{x) = f{x), (1) 

f i. in a space with periodic boundary conditions. Our work 
is to a considerable extent motivated by a goal to solve 
the 3D acoustic wave equation with position dependent 
material properties, and its related inverse problem a, a, 
to interesting accuracy in acceptable time on current (2015) 
high-end laptops. 

However, the classes used to solve this problem 
are designed with additional topologies, geometries, and 
applications in mind. These classes are Lattice, 
LatticeFunction, and LatticeOperator. A spe¬ 
cific application from Quantum Mechanics a has been 
refactored to extend these classes. 
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II. The Lattice class 

This class is intended to handle the most basic properties 
and operations of a discretized model. We divide them into 
topological and geometrical aspects of the model. The most 
basic properties of a dicrete model are the dimensionality of 
space, and how we approximate a continuous space with a 
number of sites in each direction (referred to as its shape). 
The code snippet 

LI = Lattice (shape= (2**13, )) 

L2 = Lattice (bC= ('P' , 'A')) 

L3 = Lattice (shape= (2**8, 2**8, 2**7)) 

demonstrate how three Lattice instances can be defined, 
LI with a one-dimensional lattice of 2^^ = 8192 sites, 
L2 with (by default) a two-dimensional lattice of 2^ x 2^ 
sites, and L3 with a three-dimensional lattice of 2® X 2® X 
2^ = 8 388 608 sites. In this process the instance properties 
shape, dim, and size are specified or given default values. 

A. Boundary conditions 

One additional property, bC, specifies the default boundary 
conditions in all directions. These conditions specify how 
functions defined on a finite lattice is extended beyond its 
edges, as is required when applying discrete differential 
operators or operations like Fast Fourier Transforms (FFT). 

Each specific case of bC is a property of each func¬ 
tion defined on the lattice. Hence it belongs to the class 
LatticeFunction, to be used and set by methods of 
LatticeOperator. However, since bC is often the same 
for all functions and operators in a given lattice model, 
it is convenient to provide a default property, which may 
be inherited by instances of LatticeFunction and 
LatticeOperator. 

The default value of bC is ' allP' , for periodic boundary 
conditions in all directions. Otherwise, bC must be a list 
with possible entries ' P' (for periodic extension), ' S' (for 
symmetric extension), 'A' (for antisymmetric extension), 
' F' (for extension with fixed provided values), and ' Z' 
(for extension with zero values). 

For ' S' and 'A' the symmetry point is midway between 
two lattice points. The boundary condition can be specified 
differently in different directions, and (unless periodic ' P' ) 
differently at the two edges of a given direction (in which 
case the corresponding entry in bC must be a two-component 
list). Internally bC is either stored as [ [ ' allP' , ] ], or as 

a dim-component list of two-component lists. The Lattice 
class is equipped with a method, set_bC (bC=' allP' ), 
which returns the internal representation from a variety of 
possible inputs. 

B. Subdomains and slices 

Assume that (j){n) and (j)o{n) are two arrays defined on a 
3-dimensional lattice, with sq a constant, and that we want 


to perform the operation 

(j)o{n) = (l)o{n) + so(j){n). (2) 

Python code for this operation could be the snippet 

for nx in range(phi. shape [0]): 

for ny in range(phi. shape [1]): 

for nz in range(phi. shape [2]): 
phiO[nx, ny, nz] = \ 

phiO[nx, ny, xz] + \ 
sO*phi[nx, ny, nz] 

This code is lengthy (hence error-prone) and runs slowly, 
because all for-loops are executed in native Python. The 
numpy code for the same operation is simply 

phiO -l-= sO*phi 

wherein all loop operations are delegated to numpy (and 
maybe further translated to optimized BLAS operations)Q 
The similar operation corresponding to 

(j3o{n) = (jjoin)+^so{b)(j){n-b), (3) 

b 

where 6 is a non-zero integer vector, requires more care and 
coding, since there will be values of n for which n — b falls 
outside the lattice. In such cases the expression (j){n — b) 
must be related to known values of </> by use of the boundary 
conditions. Assume a case where b = (3, 0, -2), that the 
lattice have (much) more than 3 sites in all directions, and 
that the boundary conditions is given by 
bC = [ ['P' , 'P' ] , ['S' , 'A' ] , ['A' , 'S' ] ] 

We may first treat the sites n where also n — b fall inside 
the lattice: 

phiO[3:, :,-2] += sO[3,0,-2]*phi[:-3, : , 2:] 

Here the ih'ce-notation defines a rectangular subdomain of 
the lattice. For instance, the slice [3 :, :, -2] specifies the 
intersection of (i) all planes in the x-direction except the first 
3, (ii) all planes in the y-direction, and (iii) all planes in the 
z-direction except the last 2. 

Note that array positions are counted from zero, with 
negative numbers referring to distances from the end. For 
a large lattice the above operation would cover most of the 
cases, and everything if the boundary conditions were ' Z' 
in all directions. 

In our example there are three more regions to be included: 

0 < rij; < 3, and 0 < < —2, (4a) 

3 < rix < —1, and — 2 < nz < —1, (4b) 

0 < ria; < 3, and — 2 < Uz < —1- (4c) 

The case ( |4a| ) can be handled by the code 

phiO [ : 3, : , : -2 ] -l-= \ 

sO[3,0,-2]*phi[-3 : , : ,2 : ] 

using the periodic boundary condition in the x-direction. For 
the cases (|4§ and ( |4cl l two planes in the z-direction fall 
outside the lattice on the upper side. Due to the symmetric 
' S' boundary condition at this edge of the lattice, the 

*Note that the codeline phiO = phiO + s0*phi is not equivalent 
to phiO += s0*phi. In the former a new copy of phiO is made; this 
requires more memory. 


function values on these planes are related to their values 
on the last two planes inside the lattice (counted in opposite 
order). This can be handled by the code 

phiO[3:,:,-2:] += \ 

sO[3,0,-2]*phi[:-3,:,:-3:-l] 
phiO [ : 3, : , -2 : ] -l-= \ 

sO[3,0,-2]*phi[-3:,:,:-3:-l] 

For detailed information about indexing and slicing in 
numpy, consult the Indexing section of the Numpy reference 
manual Q. However, gory details like the above are best 
handled by computers. The Lattice class provides a 
method, targetNsource (b, bC=None), which yields 
all the source and target slices required for a given vector b. 
Using this, the code snippet 

for cf, dT, dS in L3 .targetNsource (b) : 

phiOldT] += cf*phi[dS] 

replaces all operations above. Here the coefficient cf is — 1 
if an odd number of antisymmetric boundary conditions are 
employed (otherwise 4-1). 

Lattice also provides a related method, 
domain (shape, shift). This returns a slice dl 
pointing to a rectangular subdomain of the lattice, of shape 
shape, shifted from the origin by an integer vector shift. 

C. Index arrays and broadcasting 

Each site of a dim-dimensional lattice is labeled by a dim¬ 
dimensional integer index vector n. To construct an array A 
defined on all points of a 3-dimensional lattice, one could 
write a code snippet similar to the following 

defA = lambda n: \ 

numpy . exp ( -numpy . dot ( n , n ) ) 
shape — {2*’;*r8, 2**8, 2**8) 

A = numpy . zeros (shape) 

for n in numpy . ndindex (shape) : 

A[n] = defA(numpy .array (n) ) 

Although this code is brief and general with respect to 
dimensionality, it is not a good way to do it. Since the f or- 
loop will be executed in native Python, the code will run too 
slow. A better way is to define three arrays nO, nl, n2, all 
of shape (2®, 2®, 2®), once and for all. We may then replace 
the code above with the snippet 
defA = lambda nO, nl, n2 : \ 

exp (-n0*n0) *exp (-nl*nl) *exp (-n2*n2) 

A = defA(nO, nl, n2) 

All loops are now implicit, and will be executed by compiled 
numpy functions. 

Further, the memory cost of permanently storing three 
large arrays can be avoided by use of the broadcasting 
facility of numpy. Since the index array nO is constant in 
the y- and z-directions, it only contains a one-dimensional 
amount of information, stored in an array of shape (2®, 1,1). 
Likewise, nl can be stored in an array of shape (1,2®,1), 
and n2 in an array of shape (1,1,2®). All these arrays 
contain the same amount of data (2® linearely stored entries). 
But, due to their different shape they will act differently 
under f.i. algebraic operations: n0*n0 will still produce an 


array of shape (2®, 1,1), and similary nl*nl an array of 
shape (1,2®, 1). However, the addition of these two results 
produces an array of shape (2®, 2®, 1). 

Finally, adding n2*n2 generates an array of the final 
shape (2®, 2®, 2®). Hence, the cost of computing and storing 
index arrays are modest. We have chosen not to include them 
as properties, but provide a method narr ( ) which computes 
them when needed. This method returns a list n of arrays, 
[n[0] , n[l] , . . 

D. Geometric properties 

The discussion above maily concerns topological proper¬ 
ties of the lattice. For most application we also need some 
geometric properties. In general these may be implemented 
by defining a dim-dimensional vector of arrays, r(n), spec¬ 
ifying the position coordinates of all sites. These coordinates 
(which should depend monotoneously on n) could also be 
dynamical, i.e. part of the equation system to be solved. 

The wide range of possibilities indicate that several ver¬ 
sions of r{n) should be implemented, with the appropriate 
version chosen when a Lattice instance is dehned. We 
have introduced a property geometry, which specifies the 
version to be used. So far, geometry can only take the 
value ' fixedRect' , wherein rectangular regions of space, 
aligned with the lattice directions, are modelled. Such regions 
can be specihed by a dim-dimensional vector ve of edge- 
lengths, plus a vector rg specifying the position of the 
“lower left” corner of the spatial region. For a given lattice 
shape parameter, this defines a lattice cell with a vector of 
sidelengths dr, such that 

dr[d] = rE [d] /shape [d]. (5) 

The position coordinate r(n) is then defined such that its 
component in the d-direction is 

r[d] = rO[d] + dr [d] * (n [d]+1/2). (6) 


A summary of all Lattice properties, with example 
values, is as provided in the table below. 


shape 

LI. shape = (8192, ) 

L3. shape = (256, 256, 128) 

dim 

L1. dim = 1 

L 3. dim = 3 

size 

LI. size = 8192 

L3.size = 8388608 

bC 

Ll.bC = [ ['allP' ]] 

L2.bC = [ ['P' , 'P' ], ['A' , 'A' 

geometry 

LI. geometry = 'fixedRect' 

rO 

Ll.rO = (0, ) 

L3.r0 = (0, 0, 0) 

rE 

Ll.rE = (1, ) 

L3.rE = (1, 1, 1) 

dr 

LI. dr = [1.19209290e-07] 
L3.dr = [0.00390625, 


0.00390625, 0.0078125] 


HI. The LatticeFunction class 

Space does not allow us to continue with an 
equally detailed discussion of all components in the 
LatticeFunction and LatticeOperator classes. 
We will instead provide examples of uses, augmented with 
general comments. 

L = Lattice (shape= (2**15, 2**15), 
rE={18, 18), r0=(-9, -9)) 
defF = lambda r: \ 

numpy. exp {-r [0] **2/2)* \ 
numpy. exp {-r[1] ** 2/2 ) 

F = LatticeFunction (L, def_F=defF) 
to = time.time{) 

F. evalFr () 
print (L.size, 

(time.time()-t0)/L. size) 


This implementation introduces three new properties: rO, 
by default a dim-dimensional tuple with entries 0, rE, by 
default a dim-dimensional tuple with entries 1, and dr, 
calculated from equation (|^. The method rvec () returns 
a list r of arrays, [r[0], r[l],...], calculated from 
equation (|^. 


E. Lattice initialization, methods, and properties 

All currently available keyword arguments and default 
values for initialization of a Lattice instance is specihed 
by the code snippet below: 

def _ init _ (shape=(128, 128), 

bC=' allP' , geometry=' fixedRect' , 
rE=None, r0=None) : 


A value of None will invoke a default initialization process, 
following the rules discussed above. The current list of 
Lattice methods, with arguments, is as follows: 


set_bC 

domain 

targetNsource 

narr 

rvec 


(bC='allP' ) 

(shape, shift) 

(b, bC=None) 

0 

0 


Here we hrst dehne a 215 X 215 lattice model with periodic 
boundary conditions, and next a gaussian function centered 
in the middle of this lattice. The parameter rE is chosen 
large enough to make the periodic extension of this function 
smooth: It acquires a discontinuity in the hrst derivative of 
magnitude 18exp(—9^/2) « 0.5 • lO”!® or smaller (i.e., 
below double precision accuracy). 

The gaussian function is not evaluated when the instance F 
is dehned, only when we execute the method F . evalFr (). 
This method evaluates the function, and stores the result in 
the array F. values. 

The wall-clock time used to perform this computation 
on a 2013 MacBook Pro with 16 Gb of memory was 
measured to 3.26 ns per point. Note that this time is mostly 
spent multiplying double precision numbers; only 2 x 2^^ 
exponential function evaluations are performed. However, if 
the the code for def_F is changed to 

def_F = lambda r: \ 

numpy . exp (-(r[0]**2 + r[l]**2)/2) 

the execution time increases to 118 ns per point. This 
increase is partly due to the fact that the exponential function 
is now evaluated 10®° times, but also because the system now 
has to deal with two very large arrays (one for the argument 


of the exponential function, and one for final result), and is 
operating very close to the limit of available memory. A more 
detailed analysis, for lattices of various sizes (total number 
of lattice points), is shown in Fig [T] 


Evaluation time per lattice point 


10-7 


10 -* 


10 


• •• Numpy (slow) 
▼ ▼▼ Numpy (fast) 
AAA C (fast) 


■ t [sec] 

^***♦*111,111X1 

..I .I .I .I 

10 ® 10 ® 107 10 ® 10 ® 


Lattice size 


Fig. 1. Comparison of NumPy and C evaluation times for a gaussian 
defined on lattices of various sizes. The fast evaluation occur when writing 
the gaussian as exp(— a:^/2) X exp(— y^/2), the slow evaluation when 
writing it as exp[— (a:^ + y^)/2]. For these cases the wall-clock and CPU 
times are essentially the same. As can be seen, there is little to gain in 
evaluation time by writing the code in a fast, compiled language like C 
(and a lot to lose in coding time). 

A. FFT and related discrete transforms 

We may apply a discrete Fourier transformation to the 
data stored in F. values. This is done by the function 
call F . FFT (). The transformed data is stored in the array 
F . f ftvalues. The inverse transform is performed by the 
function call F . iFFT () , with the transformed data being 
stored in the array F . values (overwriting any previous 
data). 

Acctually, the method FFT () (or iFFT ()) do not 
necessarily perform a regular (multidimensional) discrete 
Fourier transform fftn (or its inverse ifftn). This is 
but one of several related discrete transforms available in 
scipy. f ftpack. Other such transforms are the discrete 
cosine transform dct (suitable for functions with symmet¬ 
ric boundary conditions on both sides), the discrete sine 
transform dst (suitable for functions with antisymmetric 
boundary conditions on both sides), the fast fourier transform 
rf ft of real data, and their inverses (idct, idst, irf ft). 
The rules are 

1) If bC[0] [0] == 'allP' the transform fftn (or 
ifftn) is used. Complex data is allowed. Otherwise, 
the data is assumed to be real, and an iterated sequence 
of transforms over all axes is executed. 

2) For directions such that bC[d] [0] == 'S' the 

transform dct (or idct) is performed. 

3) For directions such that bC[d] [0] == 'A' the 

transform dst (or idst) is performed. 

4) In all other cases the transform rfft (or irfft) is 
performed. 

Note that the bC used here is a property of 
LatticeFunction. This may be different from the cor¬ 
responding property of its lattice instance. By default they 
are equal. 

The discrete transforms above are useful because they 
allow (i) differential operators to be implemented as mul¬ 
tiplication operators on the transformed functions, and (ii) 


accurate interpolation of lattice functions outside the lattice 
sites. The latter is useful for implementation of prolongations 
in multigrid methods. To assess to which extent this is a 
practical approach, we have investigated the accuracy and 
the time requirements of these transforms. The code snippet 
below illustrate how this can be done; 

shape = (2**14, 2**14) 

L = Lattice (shape=shape, bC= { 'P' , 'P' )) 

F = LatticeFunction (L) 

F. values = numpy.random . rand (* shape) 
values = numpy . copy (myF . values) 
to = time.time(); F.FFT() 
tl = time.time(); F. iFFT () 
t2 = time.time() 

err=numpy.max(numpy.abs(F . values-values) ) 
print ( (tl-tO) /L.size, (t2-tl) /L. size, err) 

The output of this code shows that the forward transform 
takes about 68 ns per lattice point, the inverse transform about 
57 ns, and that the maximum difference between the original 
and backtransformed values is 1.7 x 10“^®. I.e., the cost of a 
one-way transform is roughly the same as 20 multiplications. 
The time per site increases by almost an order of magnitude 
for a lattice of shape = (2**14, 2**15), since this is 
close to the limit of available memory. 

We have investigated the behavior above in more detail, 
for different choices of the shape and bC parameters, 
with similar results. See Fig. The crude conclusion is 
that the transformation times grow roughly linearly with 
lattice size, with a prefactor which depends only slightly on 
transformation type and lattice dimensionality. 



Lattice size 

Fig. 2. Time used to perform a discrete lattice transformations of various 
types. Each time plotted is the sum of the forward and inverse transformation 
time. ID lattices are plotted in red, 2D lattices in magenta, and 3D lattices 
in blue. Within the range of lattice sizes allowed by available memory, the 
theoretically expected logarithmic growth of transformation time with size 
is not a very distinct feature. 


B. LatticeFunction initialization, methods, and properties 

All currently available keyword arguments for initializa¬ 
tion of a LatticeFunction instance is specified by the 
argument list below; 

def _ init _ (self, lattice, def_f=None, 

def_F=None, def_g=None, def_G=None, 
bC='allP', evalf=False, evalF=False, 
evalg=False, evalG=False) : 

















As can be inferred from the above, there are several ways 
to specify a function: (i) As a function of the index arrays, 
/(n), or as a function of the position vectors, F{r). The 
discrete transform of the function also lives on a lattice, the 
dual lattice, whose sites can be labelled by a list of index 
arrays q= [q[0], q[l],..]. We denote the geometric 
version of this lattice as reciprocal space, wherein each site 
q has a reciprocal position vector k{q). Hence, the function 
can also be specified from its discrete transformation, as the 
function (iii) g{q) or (iv) G{k). 

The current list of LatticeFunction methods is as 
follows: 


qarr () 
kvec () 
evalfn () 
evalFr () 
evalgq () 
evalGk () 
FFT 0 
iFFT () 
shift {frac) 
restrict () 

prolong () 


List index vectors for the dual lattice. 
List resiprocal position vectors. 

Compute values from def_f. 
Compute values from def_F. 
Compute fftvalues from def_g. 
Compute fftvalues from def_G. 
Discrete transformation of values. 
Inverse transformation of fftvalues. 
Return the function translated by frac. 
Return the function restricted 
to a cruder lattice. 

Return the function prolonged 
to a finer lattice. 


The current list of LatticeFunction properties is as 
follows: 


lattice 

def_f 

def_F 

def_g 

def_G 

bC 

values 

fftvalues 


Related Lattice instance. 

Possible function definition (default None). 
Possible function definition (default None). 
Possible function definition (default None). 
Possible function definition (default None). 
Boundary conditions (lattice .bC). 
Array of function values. 

Array of transformed function values. 


IV. The LatticeOperator class 

Many routines in scipy. sparse . linalg do not re¬ 
quire an explicit matrix representation of the operator under 
analysis. Only some algorithm which returns the result of 
applying the operator to a given vector is needed. Such 
algorithms can be assigned to a LinearOperator in¬ 
stance, after which it functions essentially as an explicit 
matrix representation. Such algorithms should not demand 
too much memory or computation time, but do not require 
any explicitly known sparse representation of the operator. 
F.i., any computational process involving a fixed number of 
multiplication, additions and fast fourier transformations will 
have a memory requirement which scales linearly with the 
lattice size, and a time requirement which (for large systems) 
also scales roughly linearly with lattice size. 

The LinearOperator class requires an input vector 
of shape (M, ) or (M, 1), and an output vector of shape 

(N, ). For higher-dimensional lattices this does not match 

the natural construction of lattice operators, which we do 
not want to interfere with. We have therefore implemented 
a general linOp (phiO) method, to be used as a universal 
matvec parameter for LinearOperator. The currently 
implemented code for this is 


phi = phi 0. reshape(self . lattice . shape) 
return numpy . ravel(self . varOp (phi )) 

This code assumes phiO to represent a scalar function. It 
will be extended to more general (vector, spinor, tensor,...) 
objects. The reshape and ravel operations above do not 
modify or move any data; they only change how the data is 
interpreted (the view of the data). 

The code above also call a specific method, 
varOp(phi). However, this is just a handle which 
should be assigned to the operator under analysis. The 
latter may either be an appropriate predefined method in 
the LatticeOperator class, or a method provided from 
outside. 

A. Explicit matrix representations 

It may be useful to inspect an explicit matrix represen¬ 
tation of a given operator on a small lattice. The method 
matrix (operator) provides such a representation: 

L = Lattice (shape= (4, ), rE={4,)) 

0 = LatticeOperator (L) 
laplace = 0 .matrix (0. laplace) 
print (laplace) 

The output from this code is 
[[-2 . 1 . 0 . 1 . ] 


1 . 

-2 . 

1 . 

0 

0 . 

1 . 

-2 . 

1 

1 . 

0 . 

1 . 

-2 


which is easily verified to have the correct form for a 3-stensil 
one-dimensional lattice Laplacian with periodic boundary 
conditions. We may redefine the lattice to have the ' Z' 
boundary condition: 

L = Lattice (shape= (4,), bC='Z', rE=(4,)) 
The output now becomes: 

[[- 2 . 1 . 0 . 0 . ] 


1 . 

-2 . 

1 . 

0 

0 . 

1 . 

-2 . 

1 

0 . 

0 . 

1 . 

-2 


When applied to a small two-dimensional lattice 
L = Lattice (shape= (2, 3) ,bC=' Z' , rE= (2, 3) ) 
the output for the correponding 5-stensil becomes 


[ [-4. 

2 . 

0 . 

0 . 

0 . 

0 

[ 2. 

-4 . 

2 . 

0 . 

0. 

0 . 

[ 0. 

2 . 

-4 . 

0 . 

0. 

0 . 

[ 0. 

0 . 

0 . 

-4 . 

2 . 

0 . 

[ 0. 

0 . 

0 . 

2 . 

-4 . 

2 . 

[ 0. 

0 . 

0 . 

0 . 

2 . 

-4 . 


We have found such applications of the matrix () method 
to be quite educating, and very useful for debugging pur¬ 
poses. 

The output matrix can also be used directly as input to 
all the standard (dense matrix) linear algebra routines in 
scipy. Lattice sizes up to about 10^ can be handled in 
this way, sufficient for most one-dimensional systems (and 
useful when comparing dense and iterative methods on small 
higher-dimensional systems). 

Methods for generating sparse matrix representations will 
also be implemented. 


B. Example of use 

An example illustrating the discussion above is provided 
by the code snippet; 

L = Lattice (shape= (2**8, ), 

rE=(18, ), r0=(-9, )) 
defF = lambda r: numpy . exp (-r [0] **2/2 ) 

F = LatticeFunction (L, def_F=defF, 
evalF=True) 

0 = LatticeOperator (L) 

0 . varOp = 0 . laplace 
F2values = 0 . varOp (F .values) 

In this simple case it does not matter if F2values is 
computed by use of 0 . linOp, 0 . varOp or 0 . laplace. 
The result of evaluating exp(—r^/2) can be com¬ 
pared with the exact result, (r^ — 1) exp(—r^/2). A 
good way to assess the discretization error is to compute 
max^ \A]^F{r) — AF{r)\. This is plotted in Fig. for a 
range of square lattices. 

C. The lattice Laplace operator 

We have used a simple implementation of the lattice Lapla- 
cian in the examples above. This is the common {2d + 1)- 
stensil approximation. For periodic boundary conditions the 
implemention is very simple, as indicated by the code snippet 
below: 

def laplace (self , phi): 

Lphi = numpy . zeros_like (phi) 
for d in range (self . dim) : 

Lphi += numpy . roll (phi, 1, 

axis=d) 

Lphi += numpy. roll (phi,-1, 
axis=d) 

Lphi -= 2*phi 

return Lphi/self .dr**2 

Here the roll-function rotates the entries of the phi- 
array in the d-direction by the specified amount (±1 for the 
code above). We have investigated how fast this implemen¬ 
tation is.. The results is plotted in Fig. As expected, the 



Fig. 3. The maximum absolute difference between the numerical and 
exact evaluation of the Laplace operator, divided by dF, as function of the 
linear lattice size. This shows that the error scales like dF, as expected 
for these stensils. The increase in error for large linear size is probably 
due to numerical roundoff (because dF becomes very small), the decrease 
for small linears size due to incomplete sampling of errors (too few lattice 
points to compare the functions where the error is maximum). 


evaluation times scales (essentially) linearly will lattice size, 
with a prefactor which increases with the complexity of the 
stensil. But, somewhat surprisingly, the evaluation times are 
not very different from the time to make back-and-forth fast 
fourier fourier transformations. This suggests an alternative 
approach, based on fast fourier transforms. 

The roll-process is fast, with all loop operations done 
in NumPy, but requires new memory for the rolled data. 
To avoid this we have implemented a general method, 
stensOp (phi ). The essential algorithm of this is illus¬ 
trated by the snippet below: 

for b in numpy . ndindex (stensil . shape) : 

cf, dT, dS = lattice.targetNsource (b) 

phiOldT] += cf *stensil [b]*phi[dS] 

Here stensil is a (small) dim-dimensional array defining 
the operator in question. 


Evaluation times per lattice point 



Fig. 4. The times to evaluate —AlcJ), for the (standard) {2D + l)-sensil 
approximation of the Laplace operator, are plotted for various lattice sizes 
and dimensionalites. As expected, the times increases with the complexity of 
the stensil. Somewhat suiprisingly, the times are not significantly different 
from the times to perform back-and-forth fast fourier transform (or its 
discrete analogs), c.f. Fig. 
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